//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007, 2008 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

// most of this code was generated by scripts/loopderivs.py

#include <assert.h>

#include "ephi.hpp"

/** Calculate a current loop b-field and its derivatives. */
void calcCoilBField (prec_t R, prec_t a, prec_t r,
    prec_t& Ba, prec_t& dBada, prec_t& dBadr, prec_t& d2Badadr,
    prec_t& Br, prec_t& dBrda, prec_t& dBrdr, prec_t& d2Brdadr)
{
  prec_t rmod = 1.0;
  if ( r < 0 )
    {
      r = -r;
      rmod = -1.0;
    }

  // have to treat on-axis field differently
  if ( r < 1e-23 )
    {
      // math tools not up to this right now so we cheat for the dBrdr and d2Brdadr derivs
      calcCoilBField (R, a, 1e-5, Ba, dBada, dBadr, d2Badadr, Br, dBrda, dBrdr, d2Brdadr);
      Br = dBrda = 0;
      prec_t v4 = R*R;
      prec_t v3 = a*a;
      prec_t v2 = v3 + v4;
      prec_t v1 = 1/(v2*prec_t_sqrt(v2));
      prec_t v0 = 1/(v2*v2*prec_t_sqrt(v2));
      Ba = 2*MAGNETIC_CONSTANT*PREC_PI*v1*v4;
      dBada = -6*MAGNETIC_CONSTANT*PREC_PI*a*v0*v4;
      dBadr = 0;
      d2Badadr = 0;
      return;
    }

prec_t v139 = -r;
prec_t v138 = -R;
prec_t v137 = 1/R;
prec_t v136 = 4*R;
prec_t v135 = 2*R;
prec_t v134 = 1/r;
prec_t v133 = -4*r;
prec_t v132 = -3*r;
prec_t v131 = -3*R;
prec_t v130 = -2*R;
prec_t v129 = -2*r;
prec_t v128 = (R*R);
prec_t v127 = (r*r);
prec_t v126 = (a*a);
prec_t v125 = -v126;
prec_t v124 = 2*R*r;
prec_t v123 = -v127;
prec_t v122 = v138 + v139;
prec_t v121 = -2*R*r;
prec_t v120 = (1.0/(r*r));
prec_t v119 = v129 + v130;
prec_t v118 = v133 + v136;
prec_t v117 = v131 + v132;
prec_t v116 = v129 + v135;
prec_t v115 = v123 + v125 + v128;
prec_t v114 = v126 + v127 + v128;
prec_t v113 = v121 + v126 + v127 + v128;
prec_t v112 = v124 + v126 + v127 + v128;
prec_t v111 = 1/v113;
prec_t v110 = 1/v112;
prec_t v109 = 4*R*v110;
prec_t v108 = 2*R*v110;
prec_t v107 = -2*R*v110;
prec_t v106 = (1.0/(v113*v113));
prec_t v105 = 4*R*r*v110;
prec_t v104 = (1.0/(v113*v113*v113));
prec_t v103 = (1.0/(v112*v112));
prec_t v102 = -4*R*r*v110;
prec_t v101 = (1.0/(v112*v112*prec_t_sqrt(v112)));
prec_t v100 = (1.0/prec_t_sqrt(v112));
prec_t v99 = (1.0/(v112*prec_t_sqrt(v112)));
prec_t v98 = 1 + v102;
prec_t v97 = 1/v98;
prec_t v96 = prec_t_sqrt(v105);
prec_t v95 = (1.0/prec_t_sqrt(v105));
prec_t v94 = (1.0/(v98*v98));
prec_t v93 = ellipse1(v105);
prec_t v92 = ellipse2(v105, v93);
prec_t v91 = -v93;
prec_t v90 = 4*R*r*v103*v119;
prec_t v89 = 2*R*r*v103*v119;
prec_t v88 = -2*R*r*v103*v119;
prec_t v87 = 2*a*v111*v92;
prec_t v86 = 2*r*v111*v92;
prec_t v85 = -2*a*v111*v92;
prec_t v84 = -2*r*v111*v92;
prec_t v83 = 4*a*r*v106*v92;
prec_t v82 = v109 + v90;
prec_t v81 = -4*a*r*v106*v92;
prec_t v80 = v108 + v89;
prec_t v79 = v107 + v88;
prec_t v78 = -v93*v95;
prec_t v77 = v91 + v92;
prec_t v76 = v111*v115*v92;
prec_t v75 = 2*a*v106*v116*v92;
prec_t v74 = v111*v114*v92;
prec_t v73 = -2*a*v106*v116*v92;
prec_t v72 = -2*a*v106*v115*v92;
prec_t v71 = -2*a*v106*v114*v92;
prec_t v70 = v106*v115*v116*v92;
prec_t v69 = v106*v114*v116*v92;
prec_t v68 = -2*a*v104*v115*v118*v92;
prec_t v67 = -2*a*v104*v114*v118*v92;
prec_t v66 = v92*v95*v97;
prec_t v65 = v76 + v93;
prec_t v64 = v74 + v91;
prec_t v63 = 2*a*r*v110*v111*v77;
prec_t v62 = -2*a*r*v110*v111*v77;
prec_t v61 = -a*v110*v111*v115*v77;
prec_t v60 = -a*v110*v111*v114*v77;
prec_t v59 = -a*v103*v111*v115*v119*v77;
prec_t v58 = -a*v106*v110*v115*v116*v77;
prec_t v57 = -a*v103*v111*v114*v119*v77;
prec_t v56 = -a*v106*v110*v114*v116*v77;
prec_t v55 = 2*MAGNETIC_CONSTANT*v100*v134*v64;
prec_t v54 = -2*MAGNETIC_CONSTANT*a*v65*v99;
prec_t v53 = -2*MAGNETIC_CONSTANT*v126*v134*v64*v99;
prec_t v52 = 2*MAGNETIC_CONSTANT*v122*v65*v99;
prec_t v51 = -2*MAGNETIC_CONSTANT*v100*v120*v64;
prec_t v50 = -2*MAGNETIC_CONSTANT*a*v100*v120*v64;
prec_t v49 = 2*MAGNETIC_CONSTANT*v122*v134*v64*v99;
prec_t v48 = 2*MAGNETIC_CONSTANT*a*v122*v134*v64*v99;
prec_t v47 = 2*MAGNETIC_CONSTANT*v120*v126*v64*v99;
prec_t v46 = -0.25*v112*v134*v137*v79*v93*v95;
prec_t v45 = -2*MAGNETIC_CONSTANT*a*v101*v117*v65;
prec_t v44 = v66 + v78;
prec_t v43 = 0.25*v112*v134*v137*v77*v80;
prec_t v42 = -2*MAGNETIC_CONSTANT*v101*v117*v126*v134*v64;
prec_t v41 = v82*v92*v94*v95;
prec_t v40 = -0.5*a*v111*v112*v134*v137*v77*v80;
prec_t v39 = 0.5*a*v111*v112*v134*v137*v77*v80;
prec_t v38 = 0.25*v112*v134*v137*v79*v92*v95*v97;
prec_t v37 = 0.25*v111*v112*v114*v134*v137*v77*v80;
prec_t v36 = 0.25*v111*v112*v115*v134*v137*v77*v80;
prec_t v35 = -0.5*a*v106*v112*v115*v134*v137*v77*v80;
prec_t v34 = -0.5*a*v106*v112*v114*v134*v137*v77*v80;
prec_t v33 = a*v110*v44*v96;
prec_t v32 = -a*v110*v44*v96;
prec_t v31 = 0.25*v112*v134*v137*v77*v80*v95*v97;
prec_t v30 = a*v103*v119*v44*v96;
prec_t v29 = -a*v103*v119*v44*v96;
prec_t v28 = -0.25*v112*v134*v137*v44*v80;
prec_t v27 = -0.25*a*v134*v137*v44*v80*v96;
prec_t v26 = 0.25*a*v134*v137*v44*v80*v96;
prec_t v25 = -0.25*v112*v134*v137*v44*v80*v96;
prec_t v24 = 0.25*v112*v134*v137*v44*v80*v96;
prec_t v23 = v25 + v43;
prec_t v22 = v33 + v60 + v71 + v87;
prec_t v21 = v32 + v61 + v72 + v85;
prec_t v20 = -a*v110*v111*v114*v23;
prec_t v19 = -a*v110*v111*v115*v23;
prec_t v18 = 2*MAGNETIC_CONSTANT*v100*v21;
prec_t v17 = 2*MAGNETIC_CONSTANT*a*v100*v134*v22;
prec_t v16 = 2*MAGNETIC_CONSTANT*v122*v21*v99;
prec_t v15 = -2*MAGNETIC_CONSTANT*a*v100*v120*v22;
prec_t v14 = 2*MAGNETIC_CONSTANT*a*v122*v134*v22*v99;
prec_t v13 = v25 + v37 + v69 + v86;
prec_t v12 = v24 + v36 + v70 + v84;
prec_t v11 = 2*MAGNETIC_CONSTANT*v100*v13*v134;
prec_t v10 = 2*MAGNETIC_CONSTANT*v100*v12;
prec_t v9 = 2*MAGNETIC_CONSTANT*a*v100*v13*v134;
prec_t v8 = -2*MAGNETIC_CONSTANT*a*v12*v99;
prec_t v7 = -2*MAGNETIC_CONSTANT*v126*v13*v134*v99;
prec_t v6 = v28 + v31 + v38 + v41 + v46;
prec_t v5 = a*v110*v6*v96;
prec_t v4 = -a*v110*v6*v96;
prec_t v3 = v19 + v27 + v29 + v35 + v4 + v40 + v58 + v59 + v63 + v68 + v73 + v83;
prec_t v2 = v20 + v26 + v30 + v34 + v39 + v5 + v56 + v57 + v62 + v67 + v75 + v81;
prec_t v1 = 2*MAGNETIC_CONSTANT*v100*v3;
prec_t v0 = 2*MAGNETIC_CONSTANT*a*v100*v134*v2;
Ba = 2*MAGNETIC_CONSTANT*v100*v65;
dBada = v18 + v54;
dBadr = rmod * (v10 + v52);
d2Badadr = rmod * (v1 + v16 + v45 + v8);
Br = rmod * 2*MAGNETIC_CONSTANT*a*v100*v134*v64;
dBrda = rmod * (v17 + v53 + v55);
dBrdr = (v48 + v50 + v9);
d2Brdadr = (v0 + v11 + v14 + v15 + v42 + v47 + v49 + v51 + v7);
}

#define EMOD (SPEED_OF_LIGHT*SPEED_OF_LIGHT*MAGNETIC_CONSTANT)

/** Calculate a current coil b-field and its derivatives. */
void calcCoilEField (prec_t R, prec_t a, prec_t r,
    prec_t& Ea, prec_t& dEada, prec_t& dEadr, prec_t& d2Eadadr,
    prec_t& Er, prec_t& dErda, prec_t& dErdr, prec_t& d2Erdadr)
{
  prec_t rmod = 1.0;
  if ( r < 0 )
    {
      r = -r;
      rmod = -1.0;
    }

  if ( r < 1e-23 )
    {
      // math tools not up to this right now so we cheat for the dBrdr and d2Brdadr derivs
      calcCoilEField (R, a, 1e-5, Ea, dEada, dEadr, d2Eadadr, Er, dErda, dErdr, d2Erdadr);
      Er = dErda = 0;

      prec_t v6 = (R*R);
      prec_t v5 = (a*a);
      prec_t v4 = v5 + v6;
      prec_t v3 = (1.0/(v4*prec_t_sqrt(v4)));
      prec_t v2 = (1.0/(v4*v4*prec_t_sqrt(v4)));
      prec_t v1 = EMOD*v3;
      prec_t v0 = -3*EMOD*v2*v5;
      Ea = EMOD*a*v3;
      dEada = v0 + v1;
      dEadr = 0;
      d2Eadadr = 0;
      return;
    }

  prec_t v155 = -r;
  prec_t v154 = -R;
  prec_t v153 = 1/R;
  prec_t v152 = 1/r;
  prec_t v151 = 2*R;
  prec_t v150 = -4*r;
  prec_t v149 = -4*R;
  prec_t v148 = -3*r;
  prec_t v147 = -3*R;
  prec_t v146 = -2*R;
  prec_t v145 = -2*r;
  prec_t v144 = -5*r;
  prec_t v143 = -5*R;
  prec_t v142 = (r*r);
  prec_t v141 = (R*R);
  prec_t v140 = (a*a);
  prec_t v139 = -7*R;
  prec_t v138 = -7*r;
  prec_t v137 = 2*R*r;
  prec_t v136 = R + r;
  prec_t v135 = v154 + v155;
  prec_t v134 = (1.0/(r*r));
  prec_t v133 = 1/PREC_PI;
  prec_t v132 = v138 + v139;
  prec_t v131 = v145 + v146;
  prec_t v130 = v149 + v150;
  prec_t v129 = v147 + v148;
  prec_t v128 = v143 + v144;
  prec_t v127 = v137 + v140 + v141 + v142;
  prec_t v126 = 1/v127;
  prec_t v125 = 4*R*v126;
  prec_t v124 = 8*R*v126;
  prec_t v123 = 2*R*v126;
  prec_t v122 = -4*R*v126;
  prec_t v121 = -2*R*v126;
  prec_t v120 = (1.0/(v127*v127*v127));
  prec_t v119 = (1.0/(v127*v127));
  prec_t v118 = 4*R*r*v126;
  prec_t v117 = -4*R*r*v126;
  prec_t v116 = (1.0/(v127*v127*v127*v127*prec_t_sqrt(v127)));
  prec_t v115 = (1.0/prec_t_sqrt(v127));
  prec_t v114 = (1.0/(v127*v127*v127*prec_t_sqrt(v127)));
  prec_t v113 = (1.0/(v127*prec_t_sqrt(v127)));
  prec_t v112 = (1.0/(v127*v127*prec_t_sqrt(v127)));
  prec_t v111 = 1 + v117;
  prec_t v110 = 4*R*v126*v136;
  prec_t v109 = -4*R*v126*v136;
  prec_t v108 = 4*R*r*v126*v136;
  prec_t v107 = 1/v111;
  prec_t v106 = -4*R*r*v126*v136;
  prec_t v105 = prec_t_sqrt(v118);
  prec_t v104 = (1.0/prec_t_sqrt(v118));
  prec_t v103 = ellipse1(v105*v105);
  prec_t v102 = (1.0/(v111*v111));
  prec_t v101 = (1.0/(v111*v111*v111));
  prec_t v100 = ellipse2(v105*v105);
  prec_t v99 = -v103;
  prec_t v98 = v106 + v151;
  prec_t v97 = v108 + v146;
  prec_t v96 = 8*R*r*v119*v131;
  prec_t v95 = 4*R*r*v119*v131;
  prec_t v94 = 2*R*r*v119*v131;
  prec_t v93 = -4*R*r*v119*v131;
  prec_t v92 = -2*R*r*v119*v131;
  prec_t v91 = 4*R*r*v119*v131*v136;
  prec_t v90 = -4*R*r*v119*v131*v136;
  prec_t v89 = v125 + v95;
  prec_t v88 = v124 + v96;
  prec_t v87 = v123 + v94;
  prec_t v86 = v122 + v93;
  prec_t v85 = v121 + v92;
  prec_t v84 = -8*R*a*r*v100*v119;
  prec_t v83 = 2*R*v103*v111;
  prec_t v82 = 16*a*v103*v119*v141;
  prec_t v81 = -v103*v104;
  prec_t v80 = 16*a*r*v103*v119*v141;
  prec_t v79 = -8*R*a*v100*v119*v136;
  prec_t v78 = v100 + v99;
  prec_t v77 = -v100*v98;
  prec_t v76 = v100*v97;
  prec_t v75 = -8*R*a*r*v100*v119*v136;
  prec_t v74 = 16*a*r*v103*v120*v130*v141;
  prec_t v73 = -8*R*a*r*v100*v120*v130*v136;
  prec_t v72 = v100*v104*v107;
  prec_t v71 = 2*EMOD*v100*v107*v113*v133;
  prec_t v70 = v110 + v118 + v91;
  prec_t v69 = v109 + v117 + v90;
  prec_t v68 = 2*R*v103*v86;
  prec_t v67 = -6*EMOD*v100*v107*v112*v133*v140;
  prec_t v66 = 2*EMOD*v100*v107*v112*v129*v133;
  prec_t v65 = -16*EMOD*R*v100*v102*v114*v133*v140;
  prec_t v64 = 2*EMOD*a*v100*v107*v112*v129*v133;
  prec_t v63 = -16*EMOD*R*r*v100*v102*v114*v133*v140;
  prec_t v62 = -6*EMOD*v100*v107*v114*v128*v133*v140;
  prec_t v61 = -16*EMOD*R*r*v100*v102*v116*v132*v133*v140;
  prec_t v60 = a*v126*v78*v98;
  prec_t v59 = v100*v70;
  prec_t v58 = v76 + v83;
  prec_t v57 = v77 + v83;
  prec_t v56 = -2*EMOD*v107*v112*v133*v140*v78;
  prec_t v55 = a*v119*v131*v78*v98;
  prec_t v54 = -2*EMOD*v107*v114*v128*v133*v140*v78;
  prec_t v53 = -(1.0/4.0)*v103*v104*v127*v152*v153*v85;
  prec_t v52 = v72 + v81;
  prec_t v51 = (1.0/4.0)*v127*v152*v153*v78*v87;
  prec_t v50 = v100*v102*v104*v89;
  prec_t v49 = -2*a*v126*v136*v78*v87;
  prec_t v48 = 2*EMOD*v100*v102*v113*v133*v89;
  prec_t v47 = 2*EMOD*a*v100*v102*v113*v133*v89;
  prec_t v46 = -6*EMOD*v100*v102*v112*v133*v140*v89;
  prec_t v45 = -16*EMOD*R*r*v100*v101*v114*v133*v140*v88;
  prec_t v44 = (1.0/4.0)*v100*v104*v107*v127*v152*v153*v85;
  prec_t v43 = a*v126*v69*v78;
  prec_t v42 = -(1.0/4.0)*v127*v152*v153*v78*v87*v98;
  prec_t v41 = (1.0/4.0)*v127*v152*v153*v78*v87*v97;
  prec_t v40 = (1.0/2.0)*EMOD*v107*v115*v133*v152*v153*v78*v87;
  prec_t v39 = (1.0/2.0)*EMOD*a*v107*v115*v133*v152*v153*v78*v87;
  prec_t v38 = -(3.0/2.0)*EMOD*v107*v113*v133*v140*v152*v153*v78*v87;
  prec_t v37 = -4*EMOD*v102*v112*v133*v140*v78*v87;
  prec_t v36 = -2*EMOD*v102*v112*v133*v140*v78*v89;
  prec_t v35 = -(1.0/2.0)*EMOD*a*v107*v113*v133*v152*v153*v58;
  prec_t v34 = -4*EMOD*a*v102*v112*v133*v58;
  prec_t v33 = -(1.0/2.0)*EMOD*v107*v115*v133*v134*v153*v58;
  prec_t v32 = (1.0/2.0)*EMOD*a*v107*v113*v133*v134*v153*v58;
  prec_t v31 = (1.0/2.0)*EMOD*v107*v113*v133*v135*v152*v153*v58;
  prec_t v30 = -(1.0/2.0)*EMOD*a*v107*v112*v129*v133*v152*v153*v58;
  prec_t v29 = -4*EMOD*a*v102*v114*v128*v133*v58;
  prec_t v28 = (1.0/4.0)*v104*v107*v127*v152*v153*v78*v87;
  prec_t v27 = -2*R*a*v105*v111*v126*v52;
  prec_t v26 = -(1.0/4.0)*v127*v152*v153*v52*v87;
  prec_t v25 = -2*R*a*v105*v111*v119*v131*v52;
  prec_t v24 = -4*EMOD*a*v101*v112*v133*v58*v88;
  prec_t v23 = (1.0/2.0)*EMOD*v102*v115*v133*v152*v153*v58*v89;
  prec_t v22 = -(1.0/2.0)*EMOD*a*v102*v113*v133*v152*v153*v58*v89;
  prec_t v21 = 4*R*a*v105*v126*v52*v87;
  prec_t v20 = -2*R*a*v105*v126*v52*v86;
  prec_t v19 = -(1.0/4.0)*v105*v127*v152*v153*v52*v87;
  prec_t v18 = -(1.0/2.0)*a*v105*v111*v152*v52*v87;
  prec_t v17 = (1.0/2.0)*v105*v111*v127*v152*v52*v87;
  prec_t v16 = v19 + v51;
  prec_t v15 = a*v126*v16*v98;
  prec_t v14 = v27 + v60 + v75 + v80;
  prec_t v13 = -2*EMOD*v107*v112*v133*v140*v16;
  prec_t v12 = (1.0/2.0)*EMOD*v107*v115*v133*v14*v152*v153;
  prec_t v11 = -(1.0/2.0)*EMOD*v107*v115*v133*v134*v14*v153;
  prec_t v10 = (1.0/2.0)*EMOD*v107*v113*v133*v135*v14*v152*v153;
  prec_t v9 = (1.0/2.0)*EMOD*v102*v115*v133*v14*v152*v153*v89;
  prec_t v8 = v17 + v42 + v59 + v68;
  prec_t v7 = v17 + v41 + v59 + v68;
  prec_t v6 = (1.0/2.0)*EMOD*v107*v115*v133*v152*v153*v8;
  prec_t v5 = -(1.0/2.0)*EMOD*a*v107*v113*v133*v152*v153*v7;
  prec_t v4 = -4*EMOD*a*v102*v112*v133*v7;
  prec_t v3 = v26 + v28 + v44 + v50 + v53;
  prec_t v2 = -2*R*a*v105*v111*v126*v3;
  prec_t v1 = v15 + v18 + v2 + v20 + v21 + v25 + v43 + v49 + v55 + v73 + v74 + v79 + v82 + v84;
  prec_t v0 = (1.0/2.0)*EMOD*v1*v107*v115*v133*v152*v153;
  Ea = 2*EMOD*a*v100*v107*v113*v133;
  dEada = v56 + v63 + v67 + v71;
  dEadr = rmod * (v39 + v47 + v64);
  d2Eadadr = rmod * (v13 + v36 + v37 + v38 + v40 + v45 + v46 + v48 + v54 + v61 + v62 + v65 + v66);
  Er = rmod * (1.0/2.0)*EMOD*v107*v115*v133*v152*v153*v57;
  dErda = rmod * (v12 + v34 + v35);
  dErdr = v23 + v31 + v33 + v6;
  d2Erdadr = v0 + v10 + v11 + v22 + v24 + v29 + v30 + v32 + v4 + v5 + v9;
}


ShapedCoilFactory::ShapedCoilFactory (int n) :
    width(n), height(n), hscale(1.0), vscale(1.0)
{
  grid.resize(width*height);
  set_all (true);
}

ShapedCoilFactory::ShapedCoilFactory (int w, int h, prec_t hs, prec_t vs) :
    width(w), height(h), hscale(hs), vscale(vs)
{
  grid.resize(width*height);
  set_all (true);
}

ShapedCoilFactory::~ShapedCoilFactory ()
{
}

void ShapedCoilFactory::set(int x, int y, bool v)
{
  grid[y*width + x] = v;
}

void ShapedCoilFactory::set_all (bool v)
{
  for ( int i = 0; i < width; i++ )
      for ( int j = 0; j < height; j++ )
          set(i, j, v);
}

void ShapedCoilFactory::set_circle ()
{
  set_all (false);
  hscale = 1.0;
  vscale = 1.0;
  for ( int i = 0; i < width; i++ )
    {
      for ( int j = 0; j < height; j++ )
        {
          vect3d offset;
          get(i, j, offset);
          if ( offset.length() <= 0.5 )
              set (i, j, true);
        }
    }
  // now we need to shrink ourselves as we went out of bounds a bit
  hscale = 1.0 - 1.0/(prec_t)width;
  vscale = 1.0 - 1.0/(prec_t)height;
}

void ShapedCoilFactory::print ()
{
  for ( int y = 0; y < height; y++ )
    {
      for ( int x = 0; x < width; x++ )
        {
          if ( grid[y * width + x] )
              printf ("@");
          else
              printf ("_");
        }
      printf ("\n");
    }
}

bool ShapedCoilFactory::get(int x, int y, vect3d& offset) const
{
  prec_t cw = hscale / (prec_t)width;
  offset.x = -hscale/2 + cw/2 + x * cw;
  prec_t ch = vscale / (prec_t)height;
  offset.y = -vscale/2 + ch/2 + y * ch;
  offset.z = 0;
  return grid[y*width + x];
}

int ShapedCoilFactory::count() const
{
  int ret = 0;
  for ( int i = 0; i < width; i++ )
      for ( int j = 0; j < height; j++ )
          ret = ret + (grid[j * width + i] ? 1 : 0);
  return ret;
}

bool ShapedCoilFactory::make (Statics& statics, const vect3d& pos, const vect3d& normal,
    prec_t radius, prec_t cradius, prec_t I, prec_t q_per_m) const
{
  prec_t realI = I / (prec_t)count();
  prec_t realq = q_per_m / (prec_t)count();
  vect3d n = normal.normal();
  prec_t rcr;
  if ( get_width() * get_hscale() < get_height() * get_vscale() )
    {
      rcr = cradius / (prec_t)get_width();
      rcr *= get_hscale();
    }
  else
    {
      rcr = cradius / (prec_t)get_height();
      rcr *= get_vscale();
    }
  for ( int i = 0; i < get_width(); i++ )
    {
      for ( int j = 0; j < get_height(); j++ )
        {
          vect3d offset;
          if ( !get(i, j, offset) )
              continue;
          vect3d p = pos + n * offset.y * cradius * 2;
          prec_t r = radius + offset.x * cradius * 2;
#if 0
          char b[128];
          printf ("make_coil: pos=%s, r=%g, realI=%f, count=%d\n", p.sprint(b), (double)r, (double)realI, count);
#endif
          StaticRing *pc = new StaticRing(realI, p, r, rcr, normal);
          pc->setCharge (realq);
          statics.addStaticElement (pc);
        }
    }
  return true;
}

void ShapedCoilFactory::add (coils_t& coils, prec_t pos, prec_t radius, prec_t cradius,
    prec_t I, prec_t q_per_m) const
{
  prec_t realI = I / (prec_t)count();
  prec_t realq = (q_per_m * radius * 2 * PREC_PI) / (prec_t)count();
  prec_t rcr;
  if ( get_width() * get_hscale() < get_height() * get_vscale() )
    {
      rcr = cradius / (prec_t)get_width();
      rcr *= get_hscale();
    }
  else
    {
      rcr = cradius / (prec_t)get_height();
      rcr *= get_vscale();
    }
  for ( int i = 0; i < get_width(); i++ )
    {
      for ( int j = 0; j < get_height(); j++ )
        {
          vect3d offset;
          if ( !get(i, j, offset) )
              continue;
          prec_t x = pos + offset.y * cradius * 2;
          prec_t r = radius + offset.x * cradius * 2;
          coils.push_back (coil_str());
          coils.back().x = x;
          coils.back().R = r;
          coils.back().I = realI;
          coils.back().q = realq;
        }
    }
}

